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Abstract 

We present a new third-order, semi-discrete, central method for approximating 
solutions to multi-dimensional systems of hyperbolic conservation laws, convection- 
diffusion equations, and related problems. Our method is a high-order extension of the 
recently proposed second-order, semi-discrete method in 

The method is derived independently of the specific piecewise polynomial recon- 
struction which is based on the previously computed cell-averages. We demonstrate 
our results, by focusing on the new third-order CWENO reconstruction presented in 



1 21]. The numerical results we present, show the desired accuracy, high resolution and 



robustness of our method. 
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1 Introduction 

Numerical methods for approximating solutions of hyperbolic conservation laws, 

-u{x,t) + —f{u{x,t)) = 0, (1.1) 
and of the related convection-diffusion equations. 



have attracted a lot of attention in recent years (see, e.g., [^, and the references therein). 



Here, u{x,t) is a conserved quantity, f{u) is a nonlinear convection flux and Q{u,Ux) is a 
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dissipation flux satisfying tlie (weak) parabolicity condition, ^Q{u, s) > 0, Vm, s. In tlie 
most general case u = {ui, ...,Un) is an n-vector in tlie rf-spatial variables, x = (xi, ...,Xd), 
and / and Q are vector-functions. 

In this paper, we focus on the class of central schemes, all of which can be viewed as an 
extension of the well-known Lax-Friedrichs (LxF) scheme, 0. The first-order LxF method 
enjoys the major advantage of simplicity over the upwind schemes (e.g., the Godunov scheme, 
1^): no (approximate) Riemann solvers or characteristic decompositions are involved in 
its construction, and therefore, its realization for complicated multi-dimensional systems is 
rather simple. At the same time, the LxF scheme suffers from excessive numerical dissipation, 
which causes a poor (smeared) resolution of discontinuities and rarefaction waves. 

A second-order, non-oscillatory central scheme was first introduces by Nessyahu and 
Tadmor in Since then, the Nessyahu- Tadmor (NT) scheme was further extended to 

higher orders of accuracy, (also see 0, |l9l)5 as well as to the multi-dimensional systems 
(0)> in i and 0, (also [0, 0, 0). 

The main ingredient in the construction of the NT method is a second-order, non- 
oscillatory, MUSCL-type [l^, piecewise linear interpolant (instead of the piecewise constant 
one, employed in the LxF scheme) in combination with the exact solver for the time evolu- 
tion. This approach allows to significantly improve the resolution of non-smooth solutions to 
hyperbolic conservation laws, (|1.1|) , while retaining the main advantage of the LxF scheme 
- simplicity. 

Unfortunately, applying the fully-discrete NT scheme (or its higher-order extensions) to 



the second-order convection-diffusion equations, (|L2|) , does not provide the desired resolution 
of discontinuities (see, e.g., [l^). This loss of resolution occurs due to the accumula- 

tion of excessive numerical dissipation, which is typical of fully-discrete central schemes with 
small time-steps At ~ (Ax)^ (see for details). 

To circumvent this difficulty, a second-order semi-discrete central scheme was introduced 
by Kurganov and Tadmor in |T^ . This scheme has smaller dissipation than the NT scheme, 
and unlike the fully-discrete central schemes, it can be efficiently used with time-steps as 
small as required by the CFL stability restriction. 

The basic idea in the construction of the second-order semi-discrete scheme was to use a 
more accurate information about the local speed of propagation of the discontinuities. One 
was then able to derive a non-staggered semi-discrete central method, by first integrating 
over non-equally spaced control volumes, out of which a new piecewise linear interpolant was 
reconstructed and finally projected on its cell-averages (without evolving in time). The final 
step, was first introduced in [|10], in a somewhat different context of transforming staggered 
methods into non-staggered methods. 

In this paper we extend the results of [|16| by introducing a new third- order, semi- discrete, 
central scheme. Our new scheme is derived in a general form which is independent of the 
reconstruction step, as long as the reconstructed interpolant is sufficiently accurate and non- 
oscillatory. In particular, we use the new third-order CWENO reconstruction proposed in 
pT| . This reconstruction provides a third-order accurate interpolant which is built from the 
given cell-averages such that it is non-oscillatory in the essentially non-oscillatory (ENO) 
sense (see 0, [0). This interpolant is written as a convex combination of two one-sided 
linear functions and one centered parabola. In smooth regions this convex combination 
guarantees the desired third-order accuracy. It automatically switches to a second-order. 
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one-sided, linear reconstruction in the presence of large gradients. Such weighted essentially 
non-oscillatory (WENO) reconstructions were first introduced in the upwind framework, 
|TT|, after which they were extended to the central framework, [^, |2T|, |22[. 

This paper is organized as follows. We start in §0 with a brief overview of central schemes 



for conservation laws. In particular we focus in § |2.1| on the CWENO reconstruction which 
we use as the building block for our third-order method below. 

We then proceed to construct our new third-order scheme. First, we deal with the fully- 
discrete, one-dimensional setup in §|^. This new fully-discrete scheme is sketched in equation 
( p.7|) . We only give the required details that are necessary to fulfill our final goal, namely, 
to derive the semi-discrete scheme. 

With the fully-discrete scheme, (|3.7|), we are ready to approach the semi-discrete limit 



in § |4.1| . Our new third-order, one-dimensional, semi-discrete scheme is then summarized 
in equation ( [4.4| ). This scheme is written in a general form which is independent of the 
reconstruction step and can also be combined with any appropriate ODE solver for carrying 
out the time evolution. In § ^.2| we then extend our semi-discrete scheme to multidimensional 
hyperbolic and (degenerate) parabohc problems. 

We end by presenting several numerical examples in in which we approximate solutions 
to hyperbolic conservation laws as well as to convection-diffusion equations. Our new method 
is shown to enjoy the expected high-accuracy as well as the robustness and the simplicity of 
the entire family of central schemes. 



2 Central Schemes for Conservation Laws 

We briefly overview the framework of central schemes for conservation laws. Consider the 
one-dimensional system ( |1.1| ). To approximate its solutions, we introduce a spatial scale, 
Ax, and integrate over the cell I{x) := | |^ — a;| < Ax/2}, 



Ut + 



Ax 



Ax 



0. 



(2.1) 



Here and below, u denotes the average of u over /, 



u{x,t) ■=-^ J u{^^t)d^. 

I{x) 

Introducing a time scale. At, integrating in time from t to t + At and sampling (|2.1|) at 
the cells [xj,Xj+i], we obtain 



+1/2 = m"+i/2 I ^)) - /(«(^i> t))] dr, 



(2.2) 



T=t'" 



where Xj := jAx, := nAt, m" := u{xj,t"-) and := u{xj,t"-). Assuming that at time 
t = we have computed the cell-averages of the approximate solution, {«"}, we would hke 
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to utilize ( |2.2| ) to compute the cell-averages at the next time level, t""*"^ = + At. To that 
extent, we introduce a piecewise-polynomial reconstruction, 

w(x,r)^^P,(a;)xi(x), (2.3) 
j 

where Xj{^) is the characteristic function of the cell Ij := I{xj), and Pj{x) is a polynomial 
which is reconstructed from the computed cell- averages, The degree of the polynomial 

depends on the desired order of accuracy of the method. Having such an approximation to 
t"), (|2.3|) , we can easily compute the RHS of (|2.2|) . The first term, M"_|_^y2) equals 



"i+l/2 



Pj(x)dx + J Pj+i{x)dx. 



For a sufficiently small time-step. At, the solution of ( |1.1| ) subject to the initial data (|2.3|) , 
prescribed at time t = t"^, will remain smooth at some neighborhood of the grid points Xj 
for t G [t", t"'"'"^]. Hence, the integrals on the RHS of ( p.2| ) can be approximated using a 
sufficiently accurate quadrature, which is determined by the overall desired accuracy of the 
method. The values at the intermediate times which will be required in the quadrature, 
can be predicted either by a Taylor expansion or using a Runge-Kutta method (consult 

For example, a piecewise-constant reconstruction, Pj{x) = -u", and a first-order quadra- 
ture, 

^ /(w(t))rft~At/(^"), 



will result with the staggered- LxF scheme (with A := At/ Ax denoting the mesh ratio), 

A piecewise linear reconstruction, Pj{x) = + {ux)^{x — Xj), with a second-order quadra- 
ture in time (such as the mid-point rule), results with the Nessyahu-Tadmor (NT) scheme. 
Applying nonlinear limiters on the discrete slopes, (m^;)", will prevent oscillations (for details, 

HI). 



see 



To obtain a third-order central scheme, one should use a third-order, piecewise parabolic 
reconstruction together with a more accurate quadrature in time, e.g., Simpson's quadrature 
rule (see for details). 



Remarks: 

1. robustness. In order to reconstruct a non-oscillatory interpolant, one typically is 
required to use nonlinear limiters. These limiters decrease the order of accuracy of the 
method at extrema and by that they play a stabilizing role (e.g., see |2^, 0, |3l|). 
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numerical dissipation and time step. When using fully-discrete central schemes 
to approximate solutions of convection-diffusion equations, (|1.2| ), the stability restric- 
tion enforce small time-steps, At ~ (Ax)^. That is why the numerical dissipation 
is accumulated and we do not obtain high resolution of discontinuities (see [|16| for 
details). 

This problem can be avoided by using semi-discrete schemes instead of the fully-discrete 
schemes. Such a second-order, central, semi-discrete scheme was introduced in |jl6l. In 



this paper we develop a third-order, central, semi-discrete scheme with small numerical 
dissipation, which can be efficiently used with the small time-steps required due to the 
second-order operators. 



3. upwind schemes. Sampling ( |2.1|) at the cells Ij, will result with upwind schemes. 
Here, one remains with the discontinuities along the interfaces and is bound to solve 
the Riemann problems there, or at least to approximate their solutions. In the scalar, 
one-dimensional case this can be easily accomplished, but the Riemann problem has 
no known solution in the general case of systems and/or several space dimensions. 

This is the reason for why central schemes can be considered as universal methods 
for solving hyperbolic conservation laws: Riemann solvers are not involved in their 
construction, and moreover, since ( p. 21) can be carried out componentwise, no charac- 
teristic decomposition is required. 

2.1 CWENO reconstruction 

The first one-dimensional, third-order central scheme in [^], implemented the non-oscillatory 
piecewise parabolic reconstruction proposed by Liu and Osher in ||25|. Since then, a vari- 
ety of simpler reconstructions has appeared in the literature. Among these, we would like 
to mention the Central-ENO reconstruction in and the Central- WENO (CWENO) re- 



construction in and [21|, which was extended to the two-dimensional setup in pO| and 



22 



Our new third-order semi-discrete method which we develop in §^ and §^ below, can be 
integrated with any third-order, non-oscillatory reconstruction. In our numerical simulations 
presented in we will use the method recently presented in [21], which we will now briefiy 
overview. 

In each cell Ij we reconstruct a quadratic polynomial as a convex combination of three 
polynomials Pi^{x),P^{x) and Pc{x), 

Pj{x) = wMx) + wMx) + WcPcix), (2.4) 

with positive weights Wi > 0,Vi G {c, a, l}, and J^i'^i — 1- The polynomials Pi^{x),P^{x) 
correspond to left and right one-sided linear reconstructions, respectively, while Pdx) is a 
parabola, centered around Xj. 

The linear functions, Pk{x) and Pi,{x), are uniquely determined by requiring them to 
conserve the one-sided cell averages {u^,u'j_^i and respectively) as 

PAx)=u] + ^^^^^^{x-xj), P,{x)=u] + ^^^^^{x-x,). (2.5) 
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The centered parabola, Pc{x), is chosen such as to satisfy, 

Pexact(x) = c^P^ix) + cMx) + (1 - Cl - c^)Pc{x), (2.6) 

with constants q's. Here, -PexactI^;) is the unique parabola that conserves the three cell 
averages, u'j and which is given by 

^EXACT(a;) =u] + Uj{x -Xj) + ^Uj{x -xjf. (2.7) 

The approximations to the point-values of u{xj,t^),Ux{xj,t"') and Uxx{xj,t"'), are denoted 
by -u", u'j, u'j and are given by 

^' 2Ax ' ^ Ax2 

In it was shown that every symmetric selection of the constants q's in ( p.6|) will provide 



the desired third-order accuracy. For example, by taking, Cl = Cr = 1/4, equations (p.5|)-l 
yield 

Pcix) = u] - ^iu]^, - 2u] + u]_,) + 

+ - + ' - • (2-8) 

In smooth regions, the coefficients Wi of the convex combination in (|2.4| ) are chosen to 
guarantee the maximum order of accuracy (in this particular case - order three), but in 
the presence of a discontinuity they are automatically switched to the best one-sided stencil 
(which generates the least oscillatory reconstruction). The weights are taken as 

Wi = ^^, where at = i e {c,n,L}. (2.9) 

i 

The constant e guarantees that the denominator does not vanish and is taken as e = 10^^. 
The value of p may be chosen to provide the highest accuracy in smooth areas and ensure the 



non-oscillatory nature of the solution near the discontinuities (consult [0, see also [|T9|, |2T|]). 
In the value p = 2 was empirically selected, and here we use the same p in most of the 
examples presented below. Finally, the smoothness indicators, I Si, are defined as 



1 1 J 



A direct computation then results with 

IS^ = {u]-u]^,)\ IS^ = {u]^,-u]f, 

ISc = §iu]^i-2u] + u]_,f + \{u]^,-u]^,r. (2.10) 

It is easy to see that in the presence of large gradients, this reconstruction switches to one 
of the second-order one-sided linear reconstructions, Pr or Pt • For more details we refer to 

m 
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3 The Fully-Discrete One-Dimensional Construction 

In this section we present the new third-order method in the fully-discrete framework. Since 
we are mainly interested in deriving the semi-discrete scheme, we will concentrate only on 
the details which are required for that task. The scheme we derive here, is a third-order 



extension of the fully-discrete second-order scheme presented in . 

Following [0, we would like to augment the integration over the Riemann fans by a 
more accurate information about the local speed of wave propagation. We start by assuming 
that in every cell, Ij, we have reconstructed a piecewise polynomial interpolant, Pj{x,t"'), 
from the previously computed cell averages, {u^}. Then, an upper bound on the speed of 
propagation of discontinuities at the cell boundaries, Xj+1/2, is given by 

a"+i/2= max p( — {u)\ (3.1) 
where p(^) denotes the spectral radius of a matrix A, i.e., p(A) := max|Aj(A)|, with Aj(A) 

i 

being its eigenvalues. We denote by ^'^'^ ^7+i/2 ^^^^ ^^"^ right intermediate values 

of u{x,t^) at Xj+i/2, i.e., 

^i+1/2 •= Pj+l{Xj+l/2,t'^), ^j+1/2 •= Pj{Xj+l/2,t'^), 

and by C{uJ_^-^^^,u^_^_^^2) ^ curve in phase space that connects uj_^-^^^ and via the 

Riemann fan. 

Remark: In most practical applications, these local maximal speeds can be easily evaluated. 
For example, in the genuinely nonlinear or linearly degenerate case one finds that (|3.1| ) 
reduces to 



(Z^-_j_-i^^2 • — max 



{K|(Vv=')-K|K«/.))}- (3.2) 

Given the piecewise polynomial interpolant at time t", {Pj(x,t"')}, and the local speeds 
of propagation, {a^_^_if2}j we construct the fully-discrete, central method in two steps, which 
are schematically described in Figure 13]^. First, we integrate over the control volumes. 



[•'^j-l/2,V^j~l/2,r\ ^ 5 J) [''^j~l/2,r^-^j+l/2,l\ ^ [t- 5 J; ^UU [■^j+i/2,i) •^j+l/2,rJ ^ [t- ; J) 

obtaining w^^^^, and wj^y^^ respectively. Due to the finite speed of propagation, the 
points a;^+i/2,i and x]^^^. 



!,r' 

^i+1/2,/ •= — a^j^^i2^t, X^+l/2,r •= ^j+1/2 + (^]+l/2^t, 



separate between smooth and non-smooth regions. That is, the solution of equation (|0 
subject to the piecewise polynomial initial data prescribed at time t = t", may be non-smooth 
only inside the intervals [x^_^_if2i ? ^^+i/2r\ ^ ^ [t", t""*"^). 

In the second step, we repeat the non-oscillatory reconstruction (this time on a nonuni- 
formly spaced grid) and project the obtained reconstruction on the original, uniform grid, 
ending up with the cell averages at the next time level t""*"^, {m"'*'^}- This last step does not 
involve time integration, and was introduced in the context of changing staggered methods 
into non-staggered methods in [|1^ . 
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Contact Author for Figure If Necessary 
Figure 3.1: Modified Central Differencing 



We now turn to the detailed description of this algorithm. Assume that the piecewise 
polynomial reconstruction in cell Ij at time is of the form 



Pj{x, r) = Aj + Bj{x - Xj) + -Cj{x - Xj) 



(3.3) 



Then a direct computation of the integrals over the control volumes, [x^+i/2i:^]+i/2r\ ^ 
[r, r+i] and [x^_i/2,., X [t\ yields 



w 



n+l 

i+1/2 



At 



-{Bj - Bj+i 



1 



a 



J+l/2 



AtAx 



+ 



^i+l/2 



12 



(3.4) 



and 



n+l 



+ 



(Ax) 
24 



^ ■ 2 
2 AtAx 



'j-l/2 '^j+l/2)-^j + 



12 



('^i-l/2 + '^i+1/2) + 



Ax - At (a" , 

V .7-1 



J-1/2 



^i+1/2 




«i-l/2) 



^j-l/2"'i+l/2 



+ a 



'j+l/2J 



(3.5) 



[f{u{x]^i/2,i,t))dt - /(u(a;^_i/2,r,^))] dt } , 



respectively. To complete these computations, one should approximate the flux integrals on 
the RHS of ( p.4| ) and ( p.5|) using, e.g., Simpson's quadrature as described in §0. 



At this stage, the approximate cell averages, }, realize the solution at t = t 



over a nonuniform grid, which is oversampled by twice the number of the original cells at 
t = t". To convert these nonuniform averages back into the original grid, we proceed along 
the lines of . 



First, from the cell averages, w'^^l, w'^^'^, given by ( p.4|) -( p3D , we reconstruct a third- 



order, piecewise polynomial, non-oscillatory interpolant (e.g., the CWENO interpolant de- 
scribed in §|2ri|), which we will denote by u'Jllll/2i^) '"^^"''^(2^); respectively. In fact, we do 
not need any high-order reconstruction w""''^(x) since it will be averaged out (consult Figure 



3.1 



We note in passing that even for a nonuniform grid data, the CWENO interpolant can be 
written explicitly (in the spirit of §2.1| ), but these details are irrelevant for the semi-discrete 
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scheme, which will be described in §|[ At that point, all that we need is to assume that such 
a reconstruction exists and that for all j it takes the form 



1 ^ 



'^'j+l/2i^) = + Bj+1/2{X - Xj+1/2) + -Cj+1/2{X - Xj+1/2) , 

w]^\x) 



(3.6) 



w 



in the non-smooth region, (x 



i'+i/2,r)' in the smooth region, (a:"_i/2,^, a;"+i/2. 



1)1 



respectively. Given ( p. 41) , (|3.5|) and ( p.6|) , we conclude by computing the new cell averages 
at time t""*"^ according to 



u 



n+l 



1 

Ax 



J-l/2v 



i+i/2,i 



^i+1/2 



w'^^ll^{x)dx + j wj^\x)dx+ I w'^^y^{x)dx 



~n+l 



~n+l 



^j-1/2 



j-l/2,r j + l/2,i 

Aa^_i/2i,_i/2 + [1 - A(a^_i/2 + a^+1/2)] + ^«"+i/2^j-+i/2 + 

Y~ (('^i-l/2)^^i-l/2 - ('^i+l/2)^^i+l/2) + 

^ ^ '^('^i-l/2)^A-l/2 + {(^]+l/2)^^j+l/2) ■ 



(3.7) 



Remark: The third-order reconstruction ( p.6|) is necessary in order to guarantee the over- 
all third-order accuracy, since simple averaging over [Xj_^i,Xj^i] (without reconstruction) 
reduces the order of the resulting scheme (see ||10i ) . 



4 The Semi-Discrete Scheme 

We are now ready to derive our main result, which is the new third-order, semi-discrete, 
central scheme. First, we describe our ideas in the one-dimensional framework and then we 
extend them to multidimensional problems. 

4.1 One-Dimensional Problems 

We start with the derivation of the third-order semi-discrete scheme for one-dimensional 
(systems of) hyperbolic conservation laws. Using the fully-discrete scheme obtained in §^ 
the semi-discrete approximation can be directly written as the limit 

—Uiit) = lim ^. 4.1 

Substituting ( p.7[ ) into ( |4.1| ) results with 



^ = i™o{^''^~V2Vl/2-^(«^l/2 + «"+l/2)^r' + ^«"^^^^ 

+ ^«^-^")|- (4-2) 
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In the limit as At 0, all the Riemann fans have zero widths and therefore, 



A 



■i+i/2 



w 



n+1 



A 



■i-i/2 



w 



n+1 

i-1/2- 



(4.3) 



Using (|3.3| ) we can also obtain 



Ax 



U[X 



Ax „ (Ax) 



(Ax) 



pt) PA^j+i/2,t) = A, + + 



-C.J =: u. 



i+l/2(^)- 



Finally, plugging (|3.4|) , (|3.5| ) and ( [4 .31) into ([4.2|) we compute the time limit explicitly, ending 
up with our new semi-discrete scheme. 



duj 



2Ax 



/(<1/2W) + /(Vl/2(^)) - /(<-l/2(^)) - /(Vl/2(^)), 

ai-l/2(^) 



+ 



(4.4) 



2Ax 



2Ax 



u 



Vl/2(^) - ^•+1/2^ 

with local speeds aj+i/2(t), e.g., aj+i/2(t) := max |p(^§{(uT^^/2 W)) ' 
Remarks: 

1. Our third-order scheme, ([4.4|) , admits the conservative form. 



duj 
It 



Ax 



(4.5) 



with the numerical flux 

/K+i/2(^)) + /(Vi/2(^)) ai+i/2(t) 



J+1/2 



(t) : = 



u 



i/2(^) ~ Vi/2(^) • (4-6) 



This scheme is a natural generalization of the second-order semi-discrete scheme from 
0. Moreover, the second-order scheme has exactly the same form, ( |4.5| )-( ^^ ); the 
only difference is in the more accurate computation of the intermediate values, ^^_,_i/2(^) 
and u~j^-^i^{t). It is interesting to note that also the fully-discrete, staggered, second- 
and third-order central schemes have the same structure (see 

2. Similar to the case of the second-order scheme, |l^, the non-oscillatory property of 
the piecewise parabolic reconstruction, ( |3.3| ), will guarantee the non-oscillatory nature 
of our semi-discrete scheme. But unlike the piecewise linear reconstruction utilized in 
the second-order method, a piecewise parabolic reconstruction can be only essentially 
non-oscillatory. This means that, in principle, such a reconstruction may increase the 
total variation of the computed piecewise constant solution. Our numerical examples, 
however, demonstrate that the growth of the total variation is always bounded. Such 
desirable behavior of bounded total variation in the context of central- WENO schemes, 



was already observed in |]23 
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3. We would like to stress once again the simplicity of our new method, which does not 
require any (approximate) Riemann solver or any use of the characteristic variables - 
the reconstruction of piecewise polynomial interpolant, (|3.3|), is carried out component- 
wise. In particular, unlike the standard central schemes, but similar to the second-order 
semi-discrete method in |TB[, our method is based on one grid (and not on stagger- 



ing between two grids). This can be a big advantage (compared with the traditional 
central schemes) when dealing with boundary conditions and complex geometries. 



Next, let us consider the general convection-diffusion equation, (|1 . 2|) . Similar to the 
case of the second-order semi-discrete scheme, |T6[, operator splitting is not needed. We 
can apply our third-order semi-discrete scheme, ([4.5|) - ([4.6|) , to the (degenerate) parabolic 
equation, ( |1.2D , in a straightforward manner. This results in the following scheme. 



duj _ Hj+i/2{t) - Hj^i/2{t) 
dt Ax 



(4.7) 



Here, Hj^i/2{t) is our numerical convection flux, ( [4.61 ), and Qj{t) is a high-order approxima- 
tion to the diffusion term, Q{u,Ux)^- In the examples below we use the fourth-order central 
differencing of the form 



12Ax 



where 









1 r 








■~ 12Ax l 








1 r 






n,j 


■~ UAx I 








1 r 








■~ UAx I 








1 r 






-2,j 


■~ UAx [ 



-8Q{Uj_i{t), (Mx)j_ij) + Q{Uj-2{t), Mj_2,j) 



25uj-^-2(t) — 48Mj+i(t) + 36uj — 16nj_i(t) + 3uj^2(t) 

3Mj+2(t) + 10Uj+i{t) - 18Uj + QUj^i{t) - Uj-2{t) 
Ujj^2{t) — 6'Uj_|_i(t) + 18Uj — 10Mj_i(t) — ?>Uj^2{t) 

—3uj^2(t) + 16Mj+i(t) — 36uj + 48uj_i(t) — 25uj-2(t) 



(4.8) 



(4.9) 



and {uj{t)} are point-values of the reconstructed polynomials, (3.3), i.e., Uj{t) = PAxj.t). 



4.2 Multi-Dimensional Extensions 

Without loss of generality, let us consider the two-dimensional (system of) convection- 
diffusion equations, 

Ut + f{u)x + g{u)y = Q'^iU, Ux, Uy)^ + Q^(M, M^, Uy)y, (4.10) 

where the case = = corresponds to the 2D pure hyperbolic problem. 

Suppose that we have computed an approximate solution to ( [4.10|) at some time t, and 
have reconstructed a two-dimensional piecewise polynomial, third-order, essentially non- 
oscillatory interpolant over the uniform spatial grid, {xj,yk) = {jAx, kAy). 
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Following |T6|, the 2D extension of our third-order semi-discrete scheme, ( [4.7| ),( ^l6| ), can 
be written in the following form, 



du 



'j,k 



dt 



Ax Ay 

+QUt) + Ql,{t). 



(4.11) 



Here, fc(^) -^Jfc+i/2(^) ^/-numerical convection fluxes, respectively (they 

can be viewed as a generalization of the one- dimensional flux, ([4.6|) ), 



TTX 

-"j+l/2,fc 



it) 



a 



'i+1/2, 



.kit) 



u 



l/2,kit) '^j+l/2,kit) 



(4.12) 



Hlk..l2it) 



a 



y 



it) 



u 



'tk+l/2it) ^j^fc+i/2(^) 



The numerical fluxes, (|4.12| ), are expressed in terms of the intermediate values, ^j_|_i/2 



± 



fc+i/2(^)' which are obtained from the piecewise polynomial reconstruction. The local 
speeds, CLj_^_i/2 kit) and ci^j k^i/2it), are computed, e.g., by 



a 



'j+l/2,k 



(t) := max < p 



fdf 



\Q^y^3+l/2,k 



it))).P 



-1/2,.W) 



(4.13) 



a 



y 

j-,fc+l/2 



(t) := max < p 



1/2W))'P 



=+l/2(^)) 



Finally, Q^^it) Q^j kit) are high-order, central differencing approximations to the diffusion 
terms, Q''(m, Mx,%)^ and Q^(m, Ux,%)y 

Remarks: 

1. We would like to emphasize that the problem of constructing a two-dimensional, third- 
order, non-oscillatory interpolant is highly non-trivial. Several essentially 2D recon- 
structions were proposed in |2y, |21|, Alternatively, one can use one-dimensional 
CWENO reconstruction, direction by direction, in order to compute the intermediate 
values, uf+i/2,fc(^) and uf,^^y^{t). 

Following is the recipe for the computation of w^+1/2 k (^^^ computation of other inter- 
mediate values can be carried out in the similar way). 



u 



'j+l/2,k 



Wi^P^{x j+1/2) + W^P'l{Xj+i/2) + WcPciXj+1/2), 



(4.14) 
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where the P's are the polynomials introduced in §2.1| , 



2 Ax ^ Ax2 



The weights, Wi^,w-n,Wc which are given by ( |2.9| ), and are based on the smoothness 
indicators in ( p.lO| ). 

Note that the only difference between this reconstruction and the ID reconstruction, 
( p.4|) -( pnoD , is an additional term in P^{x), —-^(uj^k+i — '2uj,k + Uj^k-i), which cor- 
responds to the second derivative in the y direction and guarantees the third-order 
accuracy of the computed intermediate values. 

This 'dimension by dimension' approach, was implemented in Example 5 below. 

2. It is straightforward to extend the two-dimensional scheme, (|4.11| ), to more space 
dimensions. In particular, the dimension-by-dimension approach is a very simple and 
promising approach for multi-dimensional problems. 



5 Numerical Examples 

We conclude the paper with a number of numerical examples. Here, in order to retain 
the overall high accuracy, the semi-discrete scheme is combined with a high-order, stable 
ODE solver to complete the spatio-temporal discretization. Numerically, we observed that 
a variety of explicit methods provide satisfactory results in the context of our semi-discrete 
scheme. 

For the inviscid problems (Examples 1, 2, 3 and 5), we used the third-order total vari- 
ation diminishing (TVD) Runge-Kutta type method introduced by Shu and Osher in |p2|| . 
However, if we apply this time-integration method or any other standard Runge-Kutta type 
method to (degenerate) parabolic problems, the time-step can be very small due to their 
strict stability restrictions. 

To overcome this difficulty, we used (in Examples 4 and 5) the third-order ODE solver 
(called DUMKA3) by Medovikov, [^. This explicit method has larger stability domains 
(compared with the standard Runge-Kutta methods), which allow larger time-steps. In 
practice, DUMKA3 works as fast as implicit methods (see |^ for details). 



We abbreviate by SD3 our third-order semi-discrete scheme, which will be combined with 
the third-order TVD Runge-Kutta type method (RK3) or with DUMKA3. 

Example 1: Linear Accuracy Test 

Consider the scalar linear hyperbolic equation 

Ut + u^ = 0, xG[0,27r], (5.1) 
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augmented with the smooth initial data, u{x, 0) = sinx, and periodic boundary conditions. 
This simple problem admits a global classical solution, which was computed at time T = 1 
with a varying number of grid points, A^. 



In Table |5.1| we check the accuracy of our third-order semi-discrete scheme, SD3, coupled 
with the RK3 ODE solver. If instead of computing the approximate convergence rate between 
two consecutive mesh refinings, one approximates the convergence rate between iV = 40 and 

= 1280, the results are 3.27 in the L^-norm and 2.91 in the L°°-norm. This clearly 
demonstrates that our scheme is third-order. 

The error is measured in terms of the pointwise values. 



M — M 



^1 := Ax \uj(T) — u{xj, T)\, \\u — u\\^ac := max \uj(T) — u{xj,T)\. 



Here, u is an approximate solution, which is realized by its values at the grid points, Xj, 

Uj{T) = Pj{xj,T), 

where the P/s are the piecewise parabolic interpolants, ( p.3|) , constructed at the final time 
t = T. 



N 


L^-error 


rate 


L°°-error 


rate 


40 


4.492e-02 




2.822e-02 




80 


1.092e-02 


2.04 


1.065e-02 


1.41 


160 


2.162e-03 


2.34 


3.426e-03 


1.64 


320 


1.811e-04 


3.58 


4.705e-04 


2.86 


640 


9.267e-06 


4.29 


2.267e-05 


4.38 


1280 


5.409e-07 


4.10 


1.171e-06 


4.27 



Table 5.1: Accuracy test for the linear advection problem, ( |5.1| ); The errors at T = 1 



Example 2: Burgers' Equation 

In this example we approximate solutions to the inviscid Burgers' equation, 

=0, xG[0,27r], (5.2) 

augmented with the smooth initial data, u{x, 0) = 0.5 -|- sinx, and periodic boundary condi- 
tions. 

The unique entropy solution of (|5.2| ) develops a shock discontinuity at the critical time 
Tc = 1. Table |5.2| shows the L^- and L°°-norms of the errors at the pre-shock time T = 0.5, 
when the solution is still smooth. Once again, when approximating the convergence rate by 
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looking at the errors for = 40 and = 1280, we get 3.25 in the L^-norm and 3.10 in 
the L°°-norm. These results indicate that our method is also third-order accurate when the 
accuracy is measured in nonlinear problems. 

In Figures ^.1| ^ ^.2| we present the approximate solutions at the post-shock time T = 2, 
when the shock is well developed. The essentially non-oscillatory nature of our scheme can 
be clearly observed. 



N 


L^-error 


rate 


L°°-error 


rate 


40 


2.370e-02 




2.225e-02 




80 


5.759e-03 


2.04 


9.053e-03 


1.30 


160 


1.161e-03 


2.31 


2.921e-03 


1.63 


320 


9.541e-05 


3.61 


3.926e-04 


2.90 


640 


4.882e-06 


4.29 


1.778e-05 


4.46 


1280 


3.044e-07 


4.00 


5.732e-07 


4.96 



Table 5.2: Accuracy test for Burgers equation, (|5.2|) ; The pre-shock errors 



Contact Author for Figure If Necessary Contact Author for Figure If Necessary 



Figure 5.1: Burgers equation, (W^ 
2, N=40. 



T = Figure 5.2: Burgers equation, ( |5.2| ); T 
2, N=80. 



Example 3: Euler Equations of Gas Dynamics 

Let us consider the one-dimensional Euler system. 



d_ 

dt 



p 


+ ^ 


m 


m 


dx 


pu^ + p 


E 


_ u{E + p) _ 



0, 



p=(7-l)-(i?- V), 



where p, u, m = pu, p and E are the density, velocity, momentum, pressure and total 
energy, respectively. We solve this system subject to Sod's Riemann initial data, proposed 



m 



u{x, 0) 



Mi= (l,0,2.5f , x<0, 
Mr = (0.125, 0,0.25f, x>0. 



The approximations to the density, velocity and pressure obtained by the SD3 scheme 
with the RK3 time discretization are presented in Figures |5.3| - |5.8| . The coefficient p in the 
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smoothness indicator, (pl9|)-( pTT0D , was taken as 0.6, which seems to be the optimal value in 
this specific example. 

We would like to stress again that our SD3 scheme does not require the characteristic 
decomposition. To improve the resolution of the contact discontinuity, which is always 
smeared while the solution to the system of Euler equations is computed by the central 
method, we implemented the Artificial Compression Method (ACM) by Harten, [§. In 
the context of central schemes, the ACM can be implemented as a corrector step to the 



component- wise approach (see for details). 



Contact Author for Figure If Necessary 

Figure 5.3: Sod problem - density. 
N=200, T=0.1644. 



Contact Author for Figure If Necessary 

Figure 5.4: Sod problem - density. 
N=400, T=0.1644. 
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Figure 5.5: Sod problem - velocity. 
N=200, T=0.1644. 
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Figure 5.6: Sod problem - velocity. 
N=400, T=0.1644. 
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Figure 5.7: Sod problem - pressure. 
N=200, T=0.1644. 



Contact Author for Figure If Necessary 

Figure 5.8: Sod problem - pressure. 
N=400, T=0.1644. 



Example 4: Convection- Diffusion Equations — the Buckley-Leverett 
Model 

In this example we solve the one- dimensional Buckley-Leverett equation, 

ut + f{u)x = e{u{u)ux)^, ev{u) > 0, (5.3) 

which can be viewed as a prototype model for the two-phase fiow in oil reservoirs. Typi- 
cally, u{u) vanishes at some values of m, and thus ( p.3|) is a degenerate parabolic equation. 
Specifically, we take 

fi^) = -^r-T-. rri, Kw)=4n(l-n), £ = 0.01, 



and consider the initial value problem with the Riemann initial data, 
u{x, 0) = 



0, 0<x<l-^, 



(5.4) 



V2 



<x <1. 



The numerical solution to this problem, obtained by the SD3 scheme augmented with the 
DUMKA3 ODE solver, is presented in Figure |5.9|. 
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Figure 5.9: Buckley-Leverett model, Figure 5.10: Buckley-Leverett model, 
(PD-(U). T=0.2. (pl)-(U), including the gravitational ef- 

fect, (O). T=0.2. 



The model, (|5.3|), becomes more complicated by adding the effects of gravitation. This 
can be obtained, e.g., by taking 



+ (1 — UY 

which is non-monotone on the interval m G [0, 1]. 

The numerical solution to this initial value problem is shown in Figure p.lO| . Note that 
the exact solution to problem (|5.3|) - (|5.4|) is not available, but our solutions seem to converge 
to the physically relevant solutions in the both cases - with gravitation or without it. 

Example 5: Incompressible Euler and Navier-Stokes equations 

In this example we consider two-dimensional viscous and inviscid incompressible flow gov- 
erned by the Navier-Stokes {u > 0) and by the Euler (z/ = 0) equations, 

ut + {u- V)u + Vp = z/A'u. (5.6) 

Here, p denotes the pressure, and u = {u, v) is the two-component divergence-free velocity 
field, satisfying 

u^ + Vy = 0. (5.7) 

In the 2D case ( ^.6|) admits an equivalent scalar formulation in terms of the vorticity, 

LOt + {uuj)^ + {vuj)y = uAu, (5.8) 

where u := Vx — Uy. The incompressibility, (|5.7|), implies that equation ( |5.8| ) can be written 
in an equivalent conservative form, 

ujt + f{uj)^ + g{uj)y = uAuj, (5.9) 

with a global convection flux, {f,g) := {ulj,vuj). A second-order, fully-discrete, staggered, 
central scheme was used to solve the two-dimensional vorticity equations in p4| . This scheme 



was proved to satisfy a maximum principle for the vorticity. (For an equivalent scheme in 
the velocity formulation, see [ll3| ). 

When applied to equation (^.9|) , our two-dimensional, third-order, semi-discrete scheme, 
(pni)-(|41^), takes the form, 

dt ~ Ax Ay (5.1Uj 
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with the numerical convection fluxes, 

Uj+l/2,k{t) 



<+lAfcW+^i+l/2,fcW 



a 



j-+l/2,fe 



it) 



UJ 



j+l/2,k 



(5.11) 



V 



'j,fc+l/2 



it) 



a 



y 

j,k+l/2 



it) 



^j,k+l/2it)+^j,k+l/2it)_ - 
^tk+l/2it) - i^lk+l/2it)\ ' 

and the local speeds, 

ai+i/2,fcW := l%+i/2,fcWI, 4fc+i/2(^) •= \^3,k+i/2it)\. (5.12) 
To approximate the linear viscosity, Atu, we used the fourth-order central differencing. 



Qj,kit) 



+ 



-^j+2,kit) + IQujj+i^kit) - 30u;j,fc(t) + 16u;j_i,fc(t) - ujj-2,kit) 

12Aa;2 

-uJj,k+2it) + 16wj,fc+i(t) - 30u;j- A:(t) + 16u;j-fc_i(t) - ujj,k-2it) 



+ 



(5.13) 



12Ay2 

To compute the intermediate values of the vorticity, we use the 'dimension by dimension' 



approach described in § |4.2| : we reconstruct the corresponding CWENO interpolants in the 
X- and ^/-directions to obtain the values of oj'^j^^i2k ^^"^ ^fk+i/2- 

Another important point in the incompressible computations is that in every time step 
one has to recover the velocities, {uj^kiVj^k}) from the known values of the vorticity, {uoj^k}- 
This can be done in many different ways (consult, e.g., and the references therein). Here 
we have used a stream-function, such that /Xip = —uj, which is obtained by solving the 
nine-points Laplacian, Aipj^k = ~^j,kit)- This provides the values of the stream-function 
with fourth-order accuracy. Its gradient, Vip, then recovers the velocity fleld. 



Uj,kit) 

'Vj,kit) 
Remarks: 



-V'j,fc+2 + ^i^j,k+l - 8lpj,k-l + 'ipj,k-2 



'j+2,k 



12Ay 



(5.14) 



12Ax 



1. Observe that in this way we retain the discrete incompressibility, namely the discrete 
velocities computed in (|5.14D satisfy 



—Uj+2,k + 8Uj+i^k — 8Uj-i^k + Uj-2,k _|_ —Vj^k+2 + 8f — SVj^k-l + Vj^k-2 _ q 



12Ax 



12Ay 
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2. The point- values of the vorticity, which are required for using the nine-points Laplacian, 



were computed from its cell averages using the 'dimension by dimension' recipe, ( 4.14 )- 

(EH- 

Finally, the intermediate values of velocities can be computed, e.g., using fourth-order 
averaging. 



•,k+\/2{t) 



16 



(5.15) 



-Vj,k+2{t) + 9t;j- fc+i(t) + 9t;j- fc-i(t) - Vj,k-2(t) 
16 



We start our numerical experiments by checking the accuracy of our scheme, ( ^.10 



( p.l5| ), augmented with the DUMKA3 time discretization. We consider the Navier-Stokes 
equations, ( ^.6|) -( pl7D with u = 0.05, subject to the smooth periodic initial data (taken from 



u{x, y,0) = — cos(x) sm{y), v{x, y, 0) = sin(x) cos(?/). 



(5.16) 



The exact solution to this problem is simply an exponential decay of the initial data, given 
by 

u{x, y,t) = — cos(x) sin(y)e~^'^*, v{x, y, t) = sin(x) cos(?/)e^^''*. 

The approximate solution with different number of grid points was computed at time 
t = 2. The errors, measured in terms of vorticity in the L°°-, L^- and L^-norms are shown 



in Table 5.3 



Nx*Ny 



32*32 
64*64 
128*128 
256*256 



L°°-error rate L^-error rate L^-error rate 



2.429e-02 - 1.791e-01 - 4.559e-02 

4.571e-03 2.41 2.814e-02 2.67 7.635e-03 2.58 

8.342e-04 2.45 3.869e-03 2.86 1.146e-03 2.74 

1.208e-04 2.79 4.966e-04 2.96 1.502e-04 2.93 



Table 5.3: Accuracy Test for the Navier-Stokes Equations. 
Errors at T = 2 



0.05. 



Next, the third-order semi-discrete scheme, ( p.l(J| )-( [5.15| ), was implemented for the peri- 
odic double shear-layer model problem taken from 0. First, we solve the Euler equations, 
( ^.6|) -( ^77D with z/ = 0, subject to the {2n, 27r)-periodic initial data, 

tanh(i(y - 7r/2), ?/ < vr, 
u{x,y,0) = v{x,y,0) = 5 ■ sm{x). (5-17) 

tanh(i(37r/2 -y), y > it, 
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Here, the "thick" shear- layer width parameter, p, is taken as ^ and the perturbation pa- 
rameter 6 = 0.05. 

The numerical results at times T = 4,6, 10 with = 64 x 64 and = 128 x 128 grid 
points are presented in Figures |5.11| - |5rT6| and ^.19H5.20| . In order to compare the quality of 
the results obtained with our new method, to previous results, we plot in Figures |5.171 - |5.l8 
the results obtained for the same double shear-layer problem with the second-order central 
scheme proposed in 
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Compared with the second-order method, the new third-order 
method, can better resolve the large gradients. Since we are using only an essentially non- 
oscillatory reconstruction, some oscillations are created with the third-order method (and 
not with the "fully" non-oscillatory second-order method). 

Finally, we solve the Navier-Stokes (N-S) equations, (|5.6|) - (|5.7|) with u = 0.01, augmented 
with the "thick" shear-layer periodic initial data, ( |5.17D . 

The numerical results at time T = 10 with iV = 64 x 64 and A^= 128 x 128 grid points 



are presented in Figures |5.21|-p?^ 



Contact Author for Figure If Necessary 

Figure 5.11: Incompressible Euler Equa- 
tions; Third-order method; T=4, 64*64 
grid. 
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Figure 5.12: Incompressible Euler Equa- 
tions; Third-order method; T=4, 128*128 
grid. 
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Figure 5.13: Incompressible Euler Equa- 
tions; Third-order method; T=6, 64*64 
grid. 
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Figure 5.14: Incompressible Euler Equa- 
tions; Third-order method; T:=6, 128*128 
grid. 
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Figure 5.15: Incompressible Euler Equa- 
tions; Third-order method; T=10, 64*64 
grid. 
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Figure 5.16: Incompressible Euler Equa- 
tions; Third-order method; T=10, 
128*128 grid. 
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Figure 5.17: Incompressible Euler Equa- 
tions; Second-order method; T=10, 64*64 
grid. 
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Figure 5.18: Incompressible Euler Equa- 
tions; Second-order method; T=10, 
128*128 grid. 
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Figure 5.19: Incompressible Euler Equa- 
tions; Third-order method; T=10, 64*64 
grid. 
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Figure 5.20: Incompressible Euler Equa- 
tions; Third-order method; T=10, 
128*128 grid. 
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Figure 5.21: Incompressible Navier-Stokes 
Equations; Third-order method; T=10, 
64*64 grid. 
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Figure 5.22: Incompressible Navier-Stokes 
Equations; Third-order method; T=10, 
128*128 grid. 
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Figure 5.23: Incompressible Navier-Stokes 
Equations; Third-order method; T=10, 
64*64 grid. 
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Figure 5.24: Incompressible Navier-Stokes 
Equations; Third-order method; T=10, 
128*128 grid. 
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